REM---- Q2inv --------- 逆写像でQ(z)のジュリア集合を描く  -------------------
REM
REM      Q(z)=z/(z^2-az+a)=q,   q=qr+i*qi は反撥的周期点, a=ar+i*ai
REM 
REM    上の方程式を解いてqの2つの逆像を求める。これらの逆像をまたqとおき，
REM    同じことをを16回繰り返して，逆写像の点( Q^(-1)(q)〜Q^(-16)(q)の全て）， 
REM                    2+4+8+ …… +2^16=2^17-2=131070（個）
REM    を全て求め，プロットする。
REM    （このような逆写像の無限個の点の集合がジュリア集合だが，現実は
REM      有限個の点だから，きれいに描けない。ファトウ集合が1点で花びら
REM      のようになっている所は特に難しい）
REM 
REM     パラメーターaの値に対応して，q は以下の3つの内どれかを選ぶ
REM        1.... 固定点 1
REM        2.... 固定点 a-1
REM        3.... 周期2の周期点  {a+1+√(a^2-2a-3)}/2
REM    (注)３つのケースは，全てが反撥的周期点のとき，3つまで選択可能
REM       （１回目白，２回目黄，3回目赤）。
REM       しかし，ジュリア集合が非有界のとき，3つ選んでも効果はほとんど
REM       ないようだ。
REM 
REM ---------------------------------- 9.21 '19 --- joe --------------------
DIM sr(17,131072),si(17,131072),co(3)
INPUT PROMPT "パラメーターaのinput (ar,ai)=":ar,ai
print"複素平面上の四角形のinput"
INPUT PROMPT"xの範囲 xl,xu=":xl,xu
INPUT PROMPT"yの範囲 yl,yu=":yl,yu
REM
READ txl,txu,tyl,tyu
DATA .1,.9,.1,.9
READ co(1),co(2),co(3)
DATA 0,6,4
REM---- 画像パラメタ設定  -----
LET wx=xu-xl
LET wy=yu-yl
LET twx=txu-txl
LET twy=tyu-tyl
LET a1=twx/wx
LET a2=twy/wy
LET b1=(xu*txl-xl*txu)/wx
LET b2=(yu*tyl-yl*tyu)/wy
SET AREA COLOR 1
PLOT AREA:txl,tyl;txu,tyl;txu,tyu;txl,tyu;txl,tyl
SET LINE COLOR 7
IF xl*xu<0 THEN PLOT LINES:b1,tyl-.02;b1,tyu+.02
IF yu*yl<0 THEN PLOT LINES:txl-.02,b2;txu+.02,b2
REM
PLOT TEXT, AT .04,.96:"Q2-Julia Set(inverse iteration),  a="
PLOT TEXT, AT .46,.96,USING"(##.#####, ##.#####)":ar,ai
PLOT TEXT, AT.04,.92:"複素平面:"
PLOT TEXT, AT .2,.92,USING"(#####.#####, #####.#####)":xl,xu
PLOT TEXT, AT .5,.92,USING"(#####.#####, #####.#####)":yl,yu
LET t=1
REM 
10 PRINT "反撥的周期点を１つ選べ"
   print"1...... 固定点１"
   print"2...... 固定点 a-1"
   print"3...... 周期2の周期点 {a+1+√(a^2-2a-3)}/2"
   INPUT PROMPT"何番ですか k=":k
   REM
   IF k=1 THEN 
      LET qr=1
      LET qi=0
      GOTO 20
   END IF 
   IF k=2 THEN 
      LET qr=ar-1
      LET qi=ai
      GOTO 20
   END IF  
   IF k=3 THEN
      LET zr=ar*ar-ai*ai-2*ar-3
      LET zi=2*ai*(ar-1)
      GOSUB 1000
      LET d1=SQR(r)*COS(s/2)
      LET d2=SQR(r)*SIN(s/2)
      LET qr=(ar+d1)/2
      LET qi=(ai+d2)/2
      GOTO 20
   ELSE
      print"再input"
      GOTO 10  
   END IF 
20 REM  -----------------------------------------------  
   PRINT"反撥的周期点 q=(";qr;",";qi;")"    
   PLOT TEXT,AT .02,.06:"反撥的周期点 q="
   IF t=3 THEN
      PLOT TEXT,AT .68,.06,USING"(##.#####, ##.#####)":qr,qi
      GOTO 30
   END IF  
   IF t=2 THEN
      PLOT TEXT,AT .44,.06,USING"(##.#####, ##.#####)":qr,qi
      GOTO 30
   END IF  
   PLOT TEXT,AT .2,.06,USING"(##.#####, ##.#####)":qr,qi 
30 REM
   REM         計算開始
   REM
   LET n=16                       !逆写像の回数
   LET ad=ar*ar+ai*ai
   LET d1=ar*ar-ai*ai
   LET i=1
   GOSUB 800
   LET sr(i,1)=ur
   LET si(i,1)=ui
   LET sr(i,2)=vr
   LET si(i,2)=vi
   FOR i=2 TO n
      FOR  j=1 TO 2^(i-1)
         LET qr=sr(i-1,j)
         LET qi=si(i-1,j) 
         GOSUB 800
         LET sr(i,2*j-1)=ur
         LET si(i,2*j-1)=ui
         LET sr(i,2*j)=vr
         LET si(i,2*j)=vi
      NEXT j
   NEXT i
    
   REM -----------
   SET POINT STYLE 1
   SET POINT COLOR co(t)
   FOR i=1 TO n
      PRINT"i=";i
      FOR j=1 TO 2^i
         LET tx=a1*sr(i,j)+b1
         IF tx<.1 OR tx>.9 THEN GOTO 40
         LET ty=a2*si(i,j)+b2
         IF ty<.1 OR ty>.9 THEN GOTO 40
         PLOT POINTS: tx,ty
40    NEXT j
   NEXT i
   REM ------------- 初期値を変える ----------------------
   INPUT PROMPT"別の初期値で計算するか？ y/n ":jo$
   IF jo$="y" OR jo$="Y" THEN
      LET t=t+1
      FOR i=1 TO n
         FOR j=1 TO 2^i
            LET sr(i,j)=0
            LET si(i,j)=0
         NEXT j
      NEXT i
      GO TO 10
   END IF 
   REM ---------------------------
   INPUT PROMPT"ファイル名，コメント書くか? y/n ":q$
   IF q$="n" OR q$="n" THEN GOTO 1100
   PRINT "ファイル名などinputせよ(全角20文字以下) (例）Q2inv-1,Q2inv-2,.."
   LINE INPUT PROMPT "ファイル名, コメントなど(全角で20文字位):":na$
   PLOT TEXT,AT.05,.02:"ファイル名など;"
   PLOT TEXT, AT .24,.02:na$
   GOTO 1100 
800 REM =================================================================
    REM  Q(z) の inverse の計算:(qr,qi) をうけて，2つのinverse images 
    REM      (ur,ui),(vr,vi)  を返す。( Q(z)=(qr,qi) の解 ）
    REM  ================================================================
    LET q1=qr*qr+qi*qi
    LET d2=qr*qr-qi*qi
    LET f1=2*(ar*qr-ai*qi)
    LET f2=2*(ar*qi+ai*qr)
    LET zr=(d1-4*ar)*d2-4*ai*qr*qi*(ar-2)+f1+1
    LET zi=2*qr*qi*(d1-4*ar)+2*ai*d2*(ar-2)+f2
    GOSUB 1000
    LET c1=SQR(r)*COS(s/2)
    LET c2=SQR(r)*SIN(s/2)
    LET ur=(ar*q1+qr*(1+c1)+qi*c2)/(2*q1)
    LET ui=(ai*q1-qi*(1+c1)+qr*c2)/(2*q1)
    LET vr=(ar*q1+qr*(1-c1)-qi*c2)/(2*q1)
    LET vi=(ai*q1-qi*(1-c1)-qr*c2)/(2*q1)
    RETURN
    REM
    REM ============================================================
1000 REM サブルーチン polor; zr,zi をinputして r,s(θのこと) を返す
     REM ========================================================-==
     LET pai=3.141592653589793
     LET zero=1.e-7
     LET r=SQR(zr*zr+zi*zi)
     IF ABS(zr)>= zero THEN GOTO 1010
     IF zi>0 THEN 1020 ELSE 1030
1020 LET s=pai/2
     RETURN
1030 IF zi<0 THEN 1040 ELSE 1050
1040 LET s=3*pai/2
     RETURN
1050 LET s=0
     RETURN
     REM ----------- 
1010 LET AT=ATN(zi/zr)
     LET s=AT+pai
     IF zr>0 THEN LET s=AT
     RETURN
     REM =====================================
1100 end